#ifndef FUNCTION
#define FUNCTION

#include <cmath>
#include <vector>

class Function {
public:
    virtual double operator() (double x) const = 0;
    virtual double derivative(double x) const{
		double h = 1e-5;
		double df = ((*this)(x + h) - (*this)(x - h)) / (2 * h);
		return df;
	}
};

std::vector<double> solveTridiagonal(
    const std::vector<double>& lower,  // 下三角元素
    const std::vector<double>& diag,   // 对角线元素
    const std::vector<double>& upper,  // 上三角元素
    const std::vector<double>& right   // 右边的常数项
) {
    int n = diag.size();
    std::vector<double> c(n), l(n), u(n);

    l[0] = diag[0];
    u[0] = upper[0] / l[0];
    c[0] = right[0] / l[0];

    for (int i = 1; i < n; ++i) {
        l[i] = diag[i] - lower[i-1] * u[i-1];
        u[i] = upper[i] / l[i];
        c[i] = (right[i] - lower[i-1] * c[i-1]) / l[i];
    }

    std::vector<double> result(n);
    result[n-1] = c[n-1];
    for (int i = n-2; i >= 0; --i) {
        result[i] = c[i] - u[i] * result[i+1];
    }

    return result;
}
#endif
